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Abstract. We investigate the evolution of the velocity distribution function of a gran- 
ular gas composed of viscoelastic particles in the homogeneous cooling state, i.e. be- 
fore clustering occurs. The deviation of the velocity distribution function from the 
Maxwellian distribution is quantified by a Sonine polynomials expansion. The first 
non- vanishing Sonine coefficient 02 (t), reveals a complex time dependence which al- 
lows to assign the granular gas the property of an age. We discuss the possibility to 
measure the age of a granular gas. 



1 Introduction 

Granular Gases as rarefied systems of granular particles in the absence of gravity 
may be exemplified by a cloud of interstellar dust. Similar as molecular gases 
Granular Gases may be described within the concepts of classical Statistical 
Mechanics, such as temperature T, velocity distribution function f{v), etc. Once 
initialized with a certain velocity distribution. Granular Gases cool down due 
to inelastic collisions of their particles. Although these systems are extremely 
simple, in principle, they reveal a variety of structure formation and much work 
has been done recently to characterize the properties of cooling Granular Gases 
(see with many references therein). 

Most of these results have been obtained under the assumption that the co- 
efficient of restitution e which characterizes the loss of energy of two colliding 
particles (see below), is a material constant. The assumption e = const., how- 
ever, does not only contradict experiments which show that e depends signifi- 
cantly on the impact velocity |2|J^ , but it contradicts even some basic mechanical 
laws 0. The simplest physically correct description of dissipative particle col- 
lisions is based on the assumption of viscoelastic material deformation during 
collisions js), which is valid for particle collisions in a certain range of impact 
velocity and is in good agreement with experimental data 

We investigate the statistical properties of Granular Gases of viscoelastic 
particles for which the dependence of the restitution coefficient on the impact 
velocity e = e(i'imp) is known Starting from a homogeneous distribution 

we study the early stage of its evolution, where no spatial structures, as clus- 
ters 1^,^ and vortexes ||l^, have emerged yet. This stage is called the homoge- 
neously cooling state. Our results show that the properties of Granular Gases 
change qualitatively if one takes into account viscoelastic material properties, 
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i.e. e — e(fimp) [0, as compared with the equivalent system, but under the 
oversimphfied assumption e ~ const.. 

In the next Section we discuss briefly the impact-velocity dependence of 
the normal restitution coefflcient. In Sec. |^ we introduce the method to study 
the velocity distribution by means of the Boltzmann equation with the Sonine 
polynomials expansion formalism, and discuss the results obtained for gases of 
particles interacting with a constant restitution coefflcient. Our main results, 
which describe the time evolution of the velocity distribution for granular gases 
of viscoelastic particles are derived in Sec. In Sec. || we discuss the concept 
of the age of a granular gas which is based on the time evolution of the velocity 
distribution. Finally, in conclusion we summarize our flndings. 



2 Two-particle interaction of viscoelastic spheres 

The microscopic dynamics of granular particles is governed by the (normal) 
restitution coefflcient e which relates the normal components of the particle 



velocities before and after a collision, 



Vi and V* 



by 



The unit vector e — r^^j/ 



gives the direction of the inter- 



center vector r^j = r,; — rj at the instant of the collision. From the conservation 
of momentum one finds the change of velocity for the colliding particles: 



(1) 



2 - ' - VV 
1 and for inelastic collisions e decreases with 



For elastic collisions one has e 
increasing degree of inelasticity. 

In literature it is frequently assumed that the restitution coefficient is a ma- 
terial constant, e — const. Experiments, e.g. [|[D, as well as theoretical investiga- 
tions 1^ show, however, that this assumption is not consistent with the nature of 
inelastic collisions, it does not agree even with a dimension analysis |4|. The im- 
pact velocity dependence of the restitution coefficient e (wimp) ^ e {^ev^^^ has 

been obtained by generalizing Hertz's contact problem to viscoelastic spheres Q. 
From the generalized Hertz equation one obtains the velocity-dependent 
restitution coefficient for viscoelastic spheres |7| 



1 - CiAa^l^ 



1/5 



2^4/5 



2/5 



with 



3/2 



(2) 



(3) 



where Y is the Young modulus, v is the Poisson ratio, and A depends on dissi- 
pative parameters of the particle material (for details see jsj). The effective mass 
and radius are defined as 



— R1R2/ {Ri 



Ri 



m^^ = 11111712/ {nil + "1-2) 



(4) 
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with i?i/2 and rriij^ being the radii and masses of the coUiding particles. The 
constants are given by 



Ci = ^3/5)^/^ ^ ^ ^ 3 2 ^ Q^9g2 . (5) 

21/552/5^(21/10) 5 ^ ^ ^ 

Equation (^) refers to the case of pure viscoelastic interaction, i.e. when the 
relative velocity \vi.je\ is not too large (to avoid plastic deformation of the par- 
ticles) and is not too small (to allow to neglect surface effects such as roughness, 
adhesion and van der Waals interactions). The dependence of e = e(|et>y |) 
(without the material dependence) was already mentioned in ||^ where heuristic 
arguments have been applied. 

In what follows we consider a granular gas of identical viscoelastic spheres of 
unit mass. 



3 Kinetics of Granular Gases: The case e = const. 

The evolution of the velocity distribution function is generally described by the 
Boltzmann-Enskog equation, which for the force-free case reads [E3UlJ,Eq|: 



^/ (t>, t) = 52(ct)(t^ y dv2 J dee (-f 12 • e) \vi2 ■ e\ 

X {xf{vr,t)f{v*2*,t) - f{vi,t)f{v2, t)} = g2{a)I{f, /) , (6) 

where tr = 2i? is the diameter of particles, 32(0') — (2 — ?7)/2(l — 77)'^ = ^ nna^ 
is packing fraction) denotes the contact value of the two-particle correlation 
function |T^, which accounts for the increasing collision frequency due to the 
excluded volume effects; e{x) is the Heaviside step-function. The velocities v" 
and ^2* refer to the precollisional velocities of the so-called inverse collision, 
which results with Vi and V2 as the after-coUisional velocities. Finally the factor 
X in the gain term appears respectively from the Jacobian of the transformation 
dvl*dv2* dvidv2 and from the relation between the lengths of the coUisional 
cylinders e \vl2 ■ e\dt = \vi2 -eldt. [jl^. For the constant restitution coefficient 
e = const, this reads x = ||l^ , [l8| , i.e., it is independent on the impact velocity 
and, therefore, independent on time. 

With the scaling Ansatz for the distribution function 

where n is the number density of the granular gas and Wo(i) is the thermal 
velocity defined by 

\nT{t)= j dv'^f{v,t),^^-nvl{t) (8) 
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the Boltzmann equation may be redueed to two independent equations: one for 
the (time- independent) scaling function /(c), and the other one for the time- 
dependence of the thermal velocity (i.e. for the temperature). Solving the equa- 
tion for the temperature one obtains T{t) = Tq/ (1 + t/r) p3[ . 

The solution of the other equation for /(c), may be found in terms of the 
Sonine polynomial expansion For the case of e = const, it reads 



/(c) = 0(c)|l+^^ap5p(c2)| , (9) 

where (f>{c) = tt"''/^ exp(— c^) is the Maxwellian distribution for the rescaled 
velocity, and the first few Sonine polynomials are 

So{x)^l, S^{x)^-x' + ^, S2ix)^^-^ + ^. (10) 

The leading (zero-order) term in Eq. (^ is the Maxwellian distribution, while 
the next-order terms, characterized by the coefficients describe the deviation 
of the distribution from the Maxwellian. For the case of e = const, the velocity 
distribution function /(c) is time-independent, therefore, the coefficients of the 
Sonine polynomials expansion are constants. 

If the inelasticity is small, one can restrict to the first non-vanishing term 
beyond the Maxwellian which has the coefficient 02 (ai — 0, according to the 
definition of temperature p4,n8|jl9]). For e = const, this coefficient reads pffl 



^ 16(l-.)(l-20 
^ 9 + 24rf + 8ed + 41e-f30(l-e)e2 ■ ^ ' 



A more accurate expression for 02 may be found in pi 



4 Kinetics of Granular Gases: Viscoelastic particles 

For viscoelastic particles the restitution coefficient e depends on the inipact ve- 
locity due to Eq. (H). Hence, the factor x in the Boltzmann equation (p|) is not 
anymore constant as for e — const, but it reads 

X = 1 + fc.Aa'^' \v,, ■ ep/^ + |cf A^a^/^ \v,, ■ e\'/' + ■■■ (12) 

Using again the scaling Ansatz (^) the rhs of the Boltzmann equation does not 
factorize into two parts, one depending only on the scaling function /(c), and 
another one, depending only on vq. Consequently, one can not obtain a set of 
decoupled equations for /(c) and for the temperature. 

Nevertheless, it is worth to substitute the generalized scaling Ansatz 

f{v,t)^-j-f{c,t) (13) 
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into the kinetic equation (||). After some algebra Eq.(P may be recast into the 
form 

f (3 + c,^)/(c,t) + B-|/(c,t)^/(/,/) (14) 
where we define the dimensionless collisional integral 

l{fj) = Jdc2jdeO{-c,2-e)\c^2-e\ [xf{cr,t)f{c;\t) - f{c^,t)f{c,,t)} (15) 
with the reduced factor x 

X = l+^-^C,S' |ci2 • e\'/' + ^^CfS'^ |ci2 • ef/' + ■■■ (16) 
which depends now on time via a quantity 

d'{t) = Ao?!'' [2r(t)]'/'" = f> [2T{t)/To]^^^° (17) 

Here S = Aa^/^[To]^/^° and Tq is the initial temperature. We also define B = 
B{t) = uo(Off2(o')<T^n, and the moments of the dimensionless collision integral 



ixp = - j dcicfli^fjj . (18) 

According to the definitions (^), (|l^) and of B, the second moment 112 defines 
the evolution of the temperature: 

^ = l92{<7)aV, I dc,c\i (/, /) = -\BT^Ji2 . (19) 

The velocity distribution function we again describe by a Sonine polynomials 
expansion as introduced in Eq. Since in contrast to the case e = const., the 
Boltzmann equation for a gas of viscoelastic particles does not factorize into a 
time dependent equation for temperature and a time-independent equation for 
the velocity distribution, the Sonine coefficients are not constants but depend 
explicitely on time, i.e. one has now 

/(c,t) = 4>{c) I 1 + f;ap(i)5p(c2) I . (20) 



Equations for ap{t) may be found by multiplying both sides of Eq.(p^) with 
cf and integrating over dci. One obtains 

00 

—p{c^)-B^^'^akVkp=^ip (21) 
k=i 

where integration by parts has been performed and where we define 

Vkp= I ^{c)cPSk{c^)dc- {cP)= f cPf{c,t)dc. (22) 
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From ( |2l| ) we see already that the granular temperature and the Sonine 
coefficients and, hence, the distribution function do not evolve independently. 

The calculation of Vkp is straightforward; the first few of them read 1^22 = 0, 
1^24 — X' "^^"^ '^'^^ moments <^c^"+^) are zero, while the even ones, (c^") may 
be expressed in terms of ak with < fc < n. Calculations show that 
implying oi = 0, according to the definition of temperature (||) (e.; 
that (c4) = ^(1 + 02). 

Now we assume that the dissipation is not large so that the deviation from 
the Maxwellian distribution may be sufficiently described by the second term in 
the expansion (|^) only, with all higher-order terms with p > 2 discarded. Then 
(Pll) is an equation for the coefficient 02. Using the above results for 1^22, J^24, 
(?) and (c'*) it is easy to show that Eq.(^l|) converts for p — 2 into identity, 
while for p = 4 it reads: 

«2 5^2(1 + 02) + 4^M4 = 0. (23) 
6 15 

With the approximation / — (j){c)[l+a2{t)S2{c^)] the time-dependent coefficients 
lip(t) may be expressed in terms of 02 due to the definition (p^). Using the 
properties of the collision integral (e.g. the conservation of the total momentum 
at collision) one can obtain relations for the /ip(i) (e.g. [Q): 

= -i J dci J dc2 J de0{-ci2 ■ e) \ci2 ■ e\ (j){ci)(f>{c2) x 
{1 + a2 [S2{cl) + S2{cl)] + al S2{cl)S2{cl)] Z\(cf + 4) 



where Aip{ci) = [ip{c*) — ip{ci)] denotes the change of some function tp{ci) in a 
direct collision. Calculations up to the second order in S yield 

2 2 

Ai2 = ^^An(5"a^ (24) 

where the coefficients A arc pure numbers: ^00 = Oj -^01 — Oj -^02 — 0, Aw ~ ujq, 
All = -^'^0, ^12 = siso'^O' "^20 = ^1, •A.21 = and ^22 = 640ooo ^1 ' '^ith 

Wo = 2v^2i/i"r (fi) Ci « 6.48562 and uji = V2^2^/^r {f) « 9.28569. 
Similarly 

2 2 

/i4 = ^^W'=a^ (25) 

k=Q n=Q 

with 600 = 0, 601 = aV2^, B02 = |\/2^, Bio = %ujo, Bii = ^uo, B12 = 

567 , , K _ 77, , R _ 149054 , , j n _ 348424 , , 

T250o'^o, O20 -j^-^i' ^21- T375o-^^i and B22 -^00000^1 

Thus, Eqs.(|23|) and (|l9|), together with Eqs.(p4[) and (g5|) form a closed set to 
find the time evolution of the temperature and of the coefficient a2. In contrast 
to the case of constant restitution coefficient where 02 — const., in a gas of 
viscoelastic particles the time evolution of temperature is coupled to the time 
evolution of the Sonine coefficient 02 . This coupling may lead to a rather peculiar 
time-evolution of the system. 
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Introducing the reduced temperature u{t) = T(t)/To we recast the set (p3|), 
into the form 

. 1 8/5 2 7 2\ 1 c 1^ 119 1547 ^\ ^ 

02 - rou^/V2 (1 + a2) + ^rou^/V4 = (27) 
where we introduce the characteristic time 

To"' = ^qoS ■ T,(0)-i = ^qoS {4V^g2iaWn^'^ (28) 

and define qo = 2i/5r(21/10)Ci/8 « 0.173318, qi = 2'^/^°{uji/ujo) ~ 1.53445, 
and = (24y/2Tr /5)qoSTo. Equation (p8|), as shown below, describes the time 
evolution of the temperature. To obtain the last equations we use the expressions 
for H2, /^4, B, and for the coefficients A. Note that the characteristic time tq is 
^ 1 times larger than the mean collision time Tc(0). 

We will find the solution to these equations as expansions in terms of the 
small dissipative parameter 5 {5'{t) = 5 ■ 2^/^°u^/^"(i)): 

u = ua + S ■ ui + S'^ ■ U2 -\ , 02 = + 5 • a2i + 5^ • 022 H (29) 



Substituting Eqs. (g9|), (|24|) and (|2^) into Eqs. (|2^), (g7|), one can solve these 
equations perturbatively for each order of S. 

Keeping only linear terms with respect to 6, one can find the analytical 



only 

solution (see fl^l for details). This reads e.g. for the coefficient 02(i): 

a2{t) = 5-021 = -^w{ty^ {Li [w{t)] - Li [w{0)]} (30) 
5 

where w{t) = exp {qo5) ^ (1 + </to)'^^ and Li(x) is the logarithmic integral. 
For t < Tq Eq.(|o|) reduces to 

02(<) = -5-/i(l-e-8*/i5-=W) (31) 

where h = 2^/1" {Bw - 5Ao) /167r = 0.415964. As it follows from Eq.(|l]), after 
a transient time of the order of few collisions per particle, i.e. for Tc(0) < i ^ tq, 
a2{t) saturates at the "steady-state" (on the time-scale ~ t'c(O)) value —h6 = 
—0.4159645. For t ^ tq it decays on a "slow" time-scale ~ tq: 

a2{t)~-5-h{tlT^)-^'^ (32) 



and the velocity distribution tends to the Maxwellian. Linear theory gives for 
the temperature for t :§> To 
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Fig. 1. Time dependence of the second coefficient of the Sonine polynomial ex- 
pansion a2{t) X 100. Time is given in units of the mean collision time Tc(0). S = 
0.1, 0.11, 0.12, . . . , 0.20 (bottom to top). 

with the constants h and qi given above. The hnear theory agrees fairly well 
with numerical solution (see Fig. |^ where the time dependence of 02 is given) 
for small S. 

For larger values of S the linear theory breaks down and we performed only 
numerical study of the equations. The results are given in Fig.^. As compared 
to the case of small S, an additional intermediate regime in the time-evolution 
of the velocity distribution is observed. The first "fast" stage of evolution takes 
place, as before, on the time scale of few collisions per particle, where maximal 
deviation from the Maxwellian distribution is achieved (Fig.||). For S > 0.15 
these maximal values of 02 are positive. Then, on the second stage (intermediate 
regime), which continues 10 — 100 collisions, 02 changes its sign and reaches a 
maximal negative deviation. Finally, on the third, slow stage, 02 (t) relaxes to 
zero on the slow time-scale ~ tq, just as for small 6. In Fig.^ (left) we show the 
first stage of the time evolution of 02 (t) for systems with large 5. At a certain 
value of the dissipative parameter S the behavior changes qualitatively, i.e. the 
system then reveals another time scale as discussed above. 
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time time 



Fig. 2. The Sonine coefficient 02 (i) for larger dissipation 5 (numerical results). 
Time is given in units of mean coUisional time Tc(0). Left: a2{t) x 100 for 5 = 
0.1,0.11,0.12, ... ,0.20 (bottom to top). Right: The plot of a-2.(t) x 100 for 5 = 0.16 
over logarithmic time shows all stages of evolution discussed in the text. 



5 The age of a granular gas 

Assume that the distribution function at time of initialization of a granular 
gas is known to be Maxwellian, e.g. from the nature of the physical process 
which gave rise to the granular gas. Then the calculation in the previous section 
describes the evolution of the velocity distribution function quantified by the 
first non-trivial Sonine coefficient 02. We have seen that this quantity evolves 
characteristically in time. Assume we know an experimental method to measure 
the velocity distribution over a certain time interval (^1,^2) and, therefore, to 
trace the time evolution of 02 in this interval. The time dependence of 02 is 
known theoretically (at least numerically) as the solution of the set of differential 
equations ( |2^p7| ), which depends on the parameters tq and S, i.e. a2 = a2{t, tq, S) 
(tq and S, in their turn, depend on the material parameters p and A). This 
suggests a method to measure the age of a granular gas, i.e., to determine the 
time of its initialization. 

Let us explain this in more detail. If we know 02 (0; ti l£ t < t2 we can 
compute the values of tq and S, which completely parameterize the dependence 
a2(t), with an accuracy, which depends on the size of the time interval ^2 — ^i- 
Hence, using the solution (cither numerical or analytical) of the set of equations 
) we can trace back the dependence 02 (i) for times t < ti. The time to when 
the curve 02 (t) cuts the abscissa corresponds to the Maxwellian distribution, i.e., 
it gives the initialization time and, therefore, the age of the gas. In Fig. || (right) 

we see that for larger dissipation there arc two such times when 02 (^0^'') = 

02 (4^^) = (0 = < tg^^ ) but only the earfier one, corresponds to the 
age (see Fig. I). Thus, we need a method to discriminate between them: If 02 (t) 
was positive for at least a part of the time interval (ti, t2), then the initialization 
time ^0 is uniquely determined (Fig. ^). If the value of 02 was negative (Fig. ^,c) 
we can trace time-backwards the dependence 02 (t) until at time <o the condition 
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02(^0) = is fulfilled. Then additional analysis is required: The value of 6 which 

was already computed determines the "steady-state" value of 02 which is reached 

after the quick relaxation in a time of the order of few mean-collision times. For 

small dissipation S it is negative and for large dissipation it is positive. Thus, if 

the steady state value of 02 is negative, the time to with 02 (fo) — corresponds 

to the initialization time, i.e. to = t'^^ Otherwise (i.e. for positive steady-state 
(2) 

02), to corresponds to tg > 0, and one has to trace a2{t) further time-backward 
in order to find the next time to which fulfills 02(^0) = to find the time of 
initialization. 

If the age of the gas has been measured according to the described method 
one can also calculate the initial temperature To, i.e. the initial energy of the gas. 
From S and Tq according to the definition of S one can estimate the combination 
of the material parameters Aa'^^^, and even the size of the particles a using 
Eq.(28), which relates tq, Tq, S and a'^g2{<y), provided the number density n 
may be measured. 

Assume an astrophysical catastrophic impact took place at a certain time 
and produced a granular gas cloud with Maxwellian velocity distribution. If one 
would be able to measure the velocity distribution function in a much later 
time interval ti < t < t2 following the described procedure one would be able 
to determine (i) the time when the impact took place, (ii) the energy of the 
impacting bodies (from the initial temperature), (iii) some material properties 
of the bodies and (iv) the grain size of the granular gas. The described analysis 
does not require the knowledge of the material properties of the particles. 



6 Conclusion 

We analyzed the time evolution of the velocity distribution function in a granular 
gas of viscoelastic particles in the homogeneous cooling state. The assumption of 
viscoelasticity is the simplest assumption for the dissipative collision of particles 
which is in agreement with mechanical laws. The collision of these particles is 
characterized by an impact-velocity dependent restitution coefficient. 

For the case of a gas of particles which interact via a constant restitution 
coefficient its evolution is completely determined by the time dependence of the 
temperature. The velocity distribution function has a simple scaling form, i.e., 
it depends only on the reduced velocity of the particles, c — v/vo{t), which is 
just the velocity measured in units of the characteristic velocity vo{t), related to 
the current temperature. The scaling form, thus, persists with time. 

Contrary, the velocity distribution function of a gas of viscoelastic particles 
does not have a simple scaling form. The deviation of the velocity distribution 
from the Maxwellian which is for the case e = const, a function of e only, i.e. time- 
independent, depends for a gas of viscoelastic particles explicitely on time, i.e. 
the velocity distribution function undergoes a time evolution. We quantify the 
deviation from the Maxwellian distribution by means of the first non-vanishing 
term of the Sonine polynomials expansion, characterized by the coefficient 02. We 
assume that inelasticity of the particles is small and, hence, higher order terms 
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Fig. 3. Illustration of the method to compute the age of the granular gas: If the velocity 
distribution and, hence, a2{t) can be measured in the interval {t\,t2), the function 
02(t) can be traced backwards in time due to the described theory. If 02(f) > for 
ti <t <t2 (Figure (a)), the condition 02 (io) = yields the age of the gas. If 02 (t) < 
for t\ < t < t2 depending on material properties there may be one time 02 (to) = 
(Figure (b)) or two times for which a2{t^Q^) = a2{t^Q^) = (Figure (c)). To discriminate 
the cases (b) and (c) one needs further consideration (see text). 



may be neglected. Contrary to the case of the constant restitution coefficient, 

where 02 = const., for a gas of viscoelastic particles 02 reveals a rather complex 
time behavior with different regimes of evolution. 

The time dependence of the distribution function, quantified by 02 (^), ex- 
hibits different stages of evolution, which allows to assign a Granular Gas the 
term "aging" . The explicite time dependence implies that the process has a def- 
inite starting point of initiahzation, i.e. earlier times do not correspond to a 
physically meaningful state of the system. At first glance it might appear cum- 
bersome that a granular gas, i.e. a cooling gas of dissipatively colliding particles, 
has a well defined initialization time. Indeed, if the temperature was the only 
characteristics of the system there would be no reason to mention some specific 
initial temperature, since the gas may exist at any temperature. In this case we 
would not be able to conclude whether the present state of a gas is the starting 
point of its evolution or an intermediate one, so that its history started long ago 
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at much higher temperature. This is the case if the particles interact via a con- 
stant restitution coefficient e = const. For gases of viscoelastic particles, however, 
the velocity distribution evolves in a way which allows to determine univocally 
the time-lag from the starting point and, hence, the age of the granular gas. 
From the knowledge of the velocity distribution in a certain time interval 
< t < ^2 we can determine not only the age of the gas but also its initial 
temperature and certain material properties of the particles. This result may be 
useful to determine the time, the energy and other system properties of catas- 
trophic impacts in astrophysical systems. 
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